/*
    Copyright 2010 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

#include<iostream>
#include "blackbody.h"
#include "wd01_Brent_PAH_grain_model.h"
#include "grain_model_hack.h"
#include "vecops.h"

//#define assert(x)

using namespace mcrx;
using namespace blitz;
using namespace std;

int main(int argc, char**argv)
{  
  thermal_equilibrium_grain g("/home/patrik/dust_data/crosssections/suv_silicate.fits",Preferences());
  //Brent_PAH_grain g("/home/patrik/dust_data/crosssections/PAH_neutral_DL07.fits");
  // we use high wavelength resolution to avoid changes due to
  // integration (log vs lin), or resampling routines. We want to
  // catch bad stuff.
  const int N=2000;
  array_1 lambda(N);
  lambda=pow(10.,(tensor::i/(1.0*N))*4-7);
  const int N2=2000;
  array_1 alambda(N2);
  alambda=pow(10.,(tensor::i/(1.0*N2))*4-7);
  g.resample(alambda,lambda);
  array_1 intens(alambda.shape());
  blackbody b(2000,.001);
  intens=b.emission(alambda);
  array_2 intens2(1,N2);
  intens2=intens(tensor::j);

  ofstream dump("grain_temp_dump.dat");
  for(int i=0;i<N;++i)
    dump << lambda(i) << "\t" << intens(i) << "\n";
  dump << "\n\n";

  //g.calculate_SED_from_intensity(intens2,mdust,dn,sed,false);
  mcrx::T_unit_map units;
  units["wavelength"]="m";
  units["length"]="m";
  units["luminosity"]="W";
  units["mass"]="Msun";
  units["time"]="yr";

  array_1 mdust(1);
  mdust=1;
  mcrx::wd01_graphite_distribution size_distribution
    ("DL07_MW3.1_60", units);
  array_1 dn(size_distribution.dn_da(g.asizes())*g.delta_size());
  array_2 sed(1,N);
  double h= g.calculate_heating(0,intens);
  double href=8.750e-24;
  cout << "\n\nTesting graphite grain emission:\n";
  cout << "Heating: " << h << ", reference value: " << href << endl;
  cout << "Difference: " << (h-href)/href << "\n\n";
  assert(abs((h-href)/href)<1e-2);

  g.calculate_SED_from_intensity(intens2,mdust,dn,sed,false);

  double lbol=integrate_quantity(sed(0,Range::all()),lambda, false);
  double lbolref=5.869e29;
  cout << "SED max index: " << maxIndex(sed(0,Range::all())) << endl;
  cout << "L_bol: " << lbol << ", reference value: " << lbolref << endl;
  cout << "Difference: " << (lbol-lbolref)/lbolref << "\n\n";
  assert(abs((lbol-lbolref)/lbolref)<1e-2);
  for(int i=0;i<N;++i)
    dump << lambda(i) << "\t" << sed(0,i) << "\n";
  dump << "\n\n";

  {
  // template emission grain
  Brent_PAH_grain pahg("/home/patrik/dust_data/crosssections/PAH_neutral_DL07.fits");
  pahg.resample(alambda,lambda);
  array_1 dn(size_distribution.dn_da(pahg.asizes())*pahg.delta_size());
  cout << "dn: " << dn << endl;
  cout << pahg.delta_size() << endl;
  cout << pahg.asizes() << endl;
  h= pahg.calculate_heating(0,intens);
  href=3.294e-26;
  cout << "\n\nTesting PAH template emission:\n";
  cout << "Heating: " << h << ", reference value: " << href << endl;
  cout << "Difference: " << (h-href)/href << endl;
  //  assert(abs((h-href)/href)<1e-2);

  pahg.calculate_SED_from_intensity(intens2,mdust,dn,sed,false);

  lbol=integrate_quantity(sed(0,Range::all()),lambda, false);
  lbolref=7.840e28;
  cout << "SED max index: " << maxIndex(sed(0,Range::all())) << endl;
  cout << "L_bol: " << lbol << ", reference value: " << lbolref << endl;
  cout << "Difference: " << (lbol-lbolref)/lbolref << "\n\n";
  //assert(abs((lbol-lbolref)/lbolref)<1e-2);
  array_1 totsig(sum(pahg.sigma()(tensor::j, tensor::i)*dn(tensor::j),tensor::j));
  cout << "total sigma is" << integrate_quantity(totsig,lambda,false) << endl;
  cout << "total heating is" << 4*constants::pi*integrate_quantity(totsig*intens,lambda,false) << endl;
  for(int i=0;i<N;++i)
    dump << lambda(i) << "\t" << sed(0,i) << '\t' << totsig(i) << "\n";
  dump << "\n\n";
  }

  //now test grain model
  Preferences prefs;
  prefs.setValue("grain_data_directory", string("/home/patrik/dust_data/crosssections"));
  prefs.setValue("wd01_parameter_set", string("MW3.1_60"));
  prefs.setValue("use_dl07_opacities", true); 
  prefs.setValue("use_cuda", false);
  prefs.setValue("n_threads", 8);
  prefs.setValue("template_pah_fraction", 0.5);
  mcrx::wd01_Brent_PAH_grain_model<mcrx::polychromatic_scatterer_policy, 
    mcrx::local_random> gm(prefs, units);
  gm.resample(alambda,lambda);
  sed=0; //sed is *added* to
  gm.calculate_SED(intens2,mdust,sed);
  lbol=integrate_quantity(sed(0,Range::all()),lambda, false);
  lbolref=5.709e30;
  cout << "\n\nTesting wd01_Brent_PAH_grain_model emission:\n";
  cout << "Model SED max index: " << maxIndex(sed(0,Range::all())) << endl;
  cout << "Model L_bol: " << lbol << ", reference value: " << lbolref << endl;
  cout << "Difference: " << (lbol-lbolref)/lbolref << endl;
  assert(abs((lbol-lbolref)/lbolref)<1e-2);
  cout << "\nTest passed.\n";
  for(int i=0;i<N;++i)
    dump << lambda(i) << "\t" << sed(0,i) << "\n";
  dump << "\n\n";
}
